Scaling, domains, and states in the four-dimensional random field Ising magnet 
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The four dimensional Gaussian random field Ising magnet is investigated numerically at zero 
temperature, using samples up to size 64 4 , to test scaling theories and to investigate the nature 
of domain walls and the thermodynamic limit. As the magnetization exponent j3 is more easily 
distinguishable from zero in four dimensions than in three dimensions, these results provide a useful 
test of conventional scaling theories. Results are presented for the critical behavior of the heat 
capacity, magnetization, and stiffness. The fractal dimensions of the domain walls at criticality 
are estimated. A notable difference from three dimensions is the structure of the spin domains: 
frozen spins of both signs percolate at a disorder magnitude less than the value at the ferromagnetic 
to paramagnetic transition. Hence, in the vicinity of the transition, there are two percolating 
clusters of opposite spins that are fixed under any boundary conditions. This structure changes 
the interpretation of the domain walls for the four dimensional case. The scaling of the effect of 
boundary conditions on the interior spin configuration is found to be consistent with the domain 
wall dimension. There is no evidence of a glassy phase: there appears to be a single transition from 
two ferromagnetic states to a single paramagnetic state, as in three dimensions. The slowing down 
of the ground state algorithm is also used to study this model and the links between combinatorial 
optimization and critical behavior. 



I. INTRODUCTION 



for the RFIM in finite dimensions. 
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As the random field Ising magnet (RFIM) is a rela- 
tively well-studied model of a disordered material, gen- 
eral questions about thermodynamic phases and transi- 
tions have been addressed using it as a model system. 
Experimental studies of random field Ising magnets are 
also available for comparison with theoretical predictions. 
The statics of the RFIM have been studied in detail 
theoretically, both analytically and numerically. It has 
been proven that there are at least two phases in di- 
mensions greater than two,El scaling arguments have-been 
constructedj3 the replica approach has been appliedj3 an 
the model has been analyzed on hierarchical lattices. 
The RFIM also has a rich numeric^ iistory, including 




extensive Monte Carlo simulatk 
perature ground state studies.Elt£ 
questions about the model remain unsettled, though, and 
the physical picture of excitations is somewhat incom- 
plete. Studying these properties of the model will be use- 
ful in building a more complete picture, especially when 
addressing questions about dynamics. 

There has been an active discussion about the na- 
ture of the phase diagram for the random field Ising 
model (RFIM). One controversy has been whether the 
transition from the ferromagnetic phase to paramagnetic 
phase, which occurs as the disorder strength or tempp.r-. 
ature is varied, is continuous in three dimensions .BuiijEj 
Recent workMO provides further strong evidence that 
the transition is second order in this case and that pre- 
viously derived scaling relations apply. However, as the 
ratio j3/u of the order parameter exponent (3 to the cor- 
relation length exponent v is very small, some scaling 
predictions are hard to verify. It is of interest to pur- 
sue this investigation in higher dimension, where f3/v is 
larger, to verify the general theoretical picture suggested 



II. SUMMARY OF RESULTS 

Numerical simulations have been carried out for the 
Gaussian RFIM on a simple hypercubic lattice in four 
dimensions. The Hamiltonian (for a review of the RFIM 
see Ref. ^ is defined over spin configurations {si = ±1}, 
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with (ij) indicating nearest neighbor sites i,j and the 
random fields hi are chosen independently from a Gaus- 
sian distribution with mean zero and variance h 2 . Here 
the energy scale is fixed by setting J = 1 in the computa- 
tions, with temperature T = 0. Exact ground states for 
this Hamiltonian ara.fauud .using a, may-flow algorithm, 
as in previous work.lBBBfflMm 

The magnetization is more useful in studying the 4D 
RFIM than the 3D RFIM, as the magnetization ex- 
ponent P is more easily distinguished from zero. The 
Binder parameter is used to locate the ferromagnetic- 
to-paramagnetic transition relatively precisely at h c — 
4.179(2). A finite-size study of the magnetization allows 
the ratio (3/v to be estimated as = 0.19(3). Besides 
its relevance to the magnetization, this ratio is impor- 
tant in studying the nature of the states and comparing 
domain wall exponents. 

The ground state energies and their dependence on 
boundary conditions can be used to study the heat ca- 
pacity and stiffness exponents of the RFIM. The stiffness 
(violation of hyperscaling) exponent is determined to be 
6 = 1.82 jjb jP.07, consistent with conventional exponent 
bounds .liijI^J Unlike the 3D case, the value for 8 is nu- 
merically distinguishable from d/2. The heat capacity 
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exponent inferred from the ground state energies is esti- 
mated as a — 0.26 ±0.05, apparently distinct from a = 
and again consistent with_the conventional disorder vari- 
ant of Widom scalingjl3£j (d — Q)v = 2 — a. 

The spatial structure of the spins for the 4D RFIM is 
found to differ from the structure for the 3D RFIM, over 
the length scales studied. In the case of three dimensions, 
the spins appear to form a nested sets of domain walls 
at criticality.E30 In 4D, the frozen spins (those invariant 
under all boundary conditions) percolate in the ferromag- 
netic phase. This implies that domain walls cannot be 
simply identified as surfaces between connected sets of 
same-sign spins. Simulations show that the frozen spins 
percolate at a value — 3.680(5). At a slightly higher 
value of h, h™, the minority spins (frozen spins of a sign 
opposite to t he magnetization) percolate. Evidence is 
given in Sec. VI that this percolation takes place even 
when the spins are coarse grained, with the critical value 
of h dependent on the scale of the coarse graining. While 
this percolation does not affect thermodynamic quanti- 
ties such as the bond part of the mean ground state 
energy, Ej, or the magnetization, the definition of the 
domain walls and the description of the spin-spin corre- 
lation function turns out to not be as straightforward as 
in the case of d = 3. 

The qualitative nature of the thermodynamic limit in 
the 4D RFIM can be addressed by studying the influ- 
ence of boundary co nditi ons on the configuration in a 
fixed window. In Sec. VII, the effect of up (+) and down 
(— ) boundary conditions at the surface of the sample are 
compared with periodic boundary conditions (P). The 
probability of the interior spins in the P configuration 
being identical to either the + or — configuration ap- 
proaches 1, as L — > oo, for all h. Taking the periodic 
boundary condition as a generic case, then, in the large 
volume limit, the interior of the ground state configura- 
tion is found in one of the two ferromagnetic states for 
h < h c or the paramagnetic state for h > h c . The prob- 
ability for the P configuration to be either + or — in 
the interior scales in a manner consistent with the 3D 
resultsES and the general case,where there are few states 
in the thermodynamic limit O 
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FIG. 1: Plot of t he sa mple average of the absolute value of 
the magnetization, \m\, of the 4D RFIM as a function of dis- 
order h and system size L, for periodic boundary conditions. 



values in which the value lies, with high confidence. In 
contrast, error bars in the figures for raw data reflect la 
statistical uncertainties computed from the standard de- 
viation, except for the Binder cumulant, where the error 
bars were computed by resampling. Plots of fitted values, 
such as estimated peak heights, include both statistical 
errors and an error bar that reflects fluctuations in values 
that result from varying the degree of the polynomial fit 
and the chosen range of the fit. 

In some of the plots, exponent values that differ from 
the "best" value from other plots are used to scale the 
data, to indicate that there is some flexibility in the ex- 
ponent values, depending on the method. All of the val- 
ues derived for the exponents from various methods are 
consistent with each other to within statistical and esti- 
mated systematic errors. Table | gives a summary of the 
numerical values of the best estimates for the exponents. 



III. MAGNETIZATION 



A. Algorithm and error bars 

The variant of the push-relabel algorithm used is the 
same as described in Rcf. Near criticality, ground 
states for samples of size 64^vere found in about 3000 s 
using 1 GHz Pentium III processors. Ground states 
for smaller samples were found using a faster, but less 
memory-efficient, version of the algorithm; using the 
same processors, near-critical samples of volume 32 4 were 
solved in approximately 60 s. 

Error bars for exponent values throughout this paper 
include both estimated systematic errors due to appar- 
ent finite size effects and errors due to statistical uncer- 
tainties; the error bars represent an estimated range of 



As the exponent (3 is more readily determined in the 
4D RFIM, compared with the 3D case, it is useful in 4D 
to study the magnetization as a first guide to the critical 
behavior and to locate the transition. The mean value of 
the absolute value of the magnetization is defined as 



\m\ = N- 1 \^2^l (2) 

i 

where the overline indicates an average over samples of 
volume N = L 4 . The magnetization is directly computed 
for each sample from the ground state with periodic 
boundary conditions. The dependence of the sample- 
averaged magnetization on disorder h is plotted in Fig. [fl. 
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TABLE I: Table of numerical estimates for the 4D Gaussian RFIM on the simple hypercubic lattice. 
Symbol Value Definition and data used 

h c 4.179 ± 0.002 Critical value of the random field for coupling J = 1. 

The critical point is determined primarily from scaling of magnetization 
distribution (e.g., Binder cumulant as shown in Fig. 

this h c is consistent with extrapolation in L of the location of peaks in the 

specific heat and the number of operations used to find the ground state 

and the value at which the probability of stiffness being zero is independent of L. 



hp 3.680 ± 0.005 Value of the random field at which the frozen spins percolate. Sec. VI. 



h™ 3.875 ± 0.005 Value of the random field at which the minority spins percolate. See Sec. V]| and Fig. 11 



P/v 0.19 ± 0.03 Ratio of magnetization exponent /3 to correlation length exponent v. 

Determined from the scaling of \m\ vs. L at criticality for 12 < L < 64. 
See Fig. |, Fig. § and Sec. 



a/v 0.31 ± 0.04 Heat capacity exponent a divided by v. 

Found from peaks C ma x(i), computed from the derivative of fit to Ej(h, L). 

See Figs. || and [j] 

(a — 1)/^ —0.94 ± 0.06 Combination of heat capacity exponent a and v. 

Found from fit to power law for the discrete estimate of dEj/d\n(L) evaluated at peak of C. 
v 0.82 ± 0.06 Correlation length exponent. 

Jointly estimated from magnetization scaling, a/v and (a — l)/v, and the scaling of the stiffness, 
with L. Consistent with the scaling of the width of the number of algorithm operations. 
P 0.16 ± 0.03 Magnetization exponent, found from j3/v and v. 
a 0.26 ± 0.05 Heat capacity exponent, found from a/v and v. 
6 1.82 ± 0.07 Violation of hyperscaling or the scaling of the stiffness at h c . . 

Found from scaling of stiffness with L and h — h c , see Sec. |v] and Fig. ^. 

d s 3.94 ± 0.06 Fractal dimension of connected domain wall at h = h c . 

Note that the result is indistinguishable from d — 4. 
See Fig. fio| and Sec. |y]. 



3.20 ± 0.12 Incongruent fractal dimension of domain wall at criticality. 

Box counting of incongruent volumes (disconnected wall). See Sec. 
Consistent with scaling of state overlap probabilities shown in Fig. 



13| and Fig, p^. 



2.94 ± 0.12 Energy "fractal dimension" at h — h c . Found from the exchange part, Ej, of the stiffness. 
See Fig. [To] and Sec. |vj 



A. Binder cumulant and h c 

One method for determining the value of h c is to 
use the Binder cumulant. The value of the cumulant, 
g = (3 — ?7i 4 /m 2 )/2, should be g = 1 in the ferromagnetic 
phase and should take on the value g = in the param- 
agnetic phase. The fixed point h c is found by the inter- 
section of the g(h) curves for various L. Some caution 
should be used with this method, as the magnetization 
exponent is small, so that the sample distribution of m is 
bimodal for even large samples near the transition. The 
assumptions of Gaussian behavior about the mean in the 
paramagnetic phase are difficult to achieve. Nonetheless, 
the plots of g{h) show a consistent behavior that indicates 
the finite-size trends in the data for the magnetization. 
The plot of g(h) for L = 4 through L = 64 is shown in 
Fig. ||. For smaller system sizes, about 5 x 10 4 ground 
states were found; for the L = 64 systems near /i c , about 



5 x 10 3 ground states were computed. The apparent in- 
tersection point is g c w 0.975, but this quantity likely 
has not converged to its scaling value. The location of 
the transition can be assigned with more certainty to the 
range h c — 4.179 ± 0.002. This is an acceptable value 
for the set of lengths used here and likely will continue 
to hold for scales that are somewhat larger. It is also 
quite consistent with the scaling of the stiffness, location 
of the specific heat peaks, and the algorithmic slowing 
down discussed in the other sections of this paper. 



B. Magnetization 

The ratio of the magnetization exponent [3 to the cor- 
relation exponent length v is computed from the effective 
finite size exponent for the magnetization. Given stan- 
dard finite size scaling, the magnetization at the transi- 
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FIG. 2: Plot of the Binder cumulant (3 - m 4 /m 2 )/2 as a 
function of h for various L. The curves are smoothed spline 
fits to indicate the trends. This plot is used to determine the 
location of the transition, h c = 4.179(2). 
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FIG. 3: Plot of estimates for (3/v found from the discrete 
derivative (/3/V)i D c = ln(|m|(/i, L 2 )/\m\(h, Li))/ln(L 2 /Li). 
The apparent convergence to a uniform value for h ~ 4.18 
implies that the value fi/v — 0.19(3) accurately describes the 
effective critical behavior for L ~ 12 — > 64. 



tion will scale as m(h c ) ~ L The local exponent /3/V 
is found by the discrete derivative of ln(|m|) with respect 
to ln(X). The results of this computation are plotted in 
Fig. y. This evaluation gives a location for the transi- 
tion that is consistent, but slightly less precise, than the 
Binder cumulant analysis. The value of the local expo- 
nent that is most consistent with a constant value gives 
the estimate 



P/v = 0.19 ±0.03. 



(3) 



In Fig. ^, the scaled magnetization (for the same samples 
used to compute g(h)) is plotted as a function of scaled 
distance to the transition, in agreement with the finite- 
size scaling form 



= L 



fm[(h - h c )L 



(4) 



where the value v = 0.83 gives the best scaling collapse, 
with fixed j3/v = 0.19 and h c = 4.179. The value of v 
in Table Q indicates the range of values found by distinct 
estimates; there is no clear best measurement of v in the 
data. 



C. Fluctuations in the magnetization 

In addition to the scaling of the mean magnetization, 
the fluctuations in the magnetization can be checked 
for consistency with finite-size scaling. The sample- 
to-samplc fluctuations Am in the total magnetization 
\M\ = iV|m| = X)i s « can be estimated, using the num- 
ber of independent volumes and the fluctuations in the 
magnetization of such volumes. Defining the finite size 
scaling variable x = [h — h^)L x l v ', at large \x\, the rele- 
vant volumes are of size £ ~ (h — h c )~ u , while at small 
larl, the volume is finite-size limited. The fluctuations in 
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FIG. 4: Scaled magnetization \m\L l3 ^ 1 ' as a function of scaled 
disorder (h - h c )L 1/v for p/v = 0.19 and v = 0.83. 



the magnetization over a volume I are M e - e d -M v and 
the number of such volumes is ng ~ (L/£) d . Taking fluc- 
tuations over each volume to be independent, one can 
write a version of the scaling as 



Mty/nt ~ L d l\h - h c ) d ^ 2 -Pf A (x), (5) 



where in the limit of large |x|, /a approaches a constant 
whose value depends on the sign of \x\. For small values 
of \x\, /a (%) ~ \x\~P +dv t 2 . This scaling form is verified 
by the data displayed in Fig. S. The data at small \x\ 
is roughly consistent with the range of power laws — (3 + 
dv/2 w 1.48 ± 0.15 plotted in Fig. | (this comparison is 
rather sensitive to the location of h c and the value of 
v.) At large \x\, A M L~ d / 2 (h - h c )> 3 ~ dl/ / 2 approaches a 
constant. 
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FIG. 5: Plot of the scaled fluctuations in the magnitude of 
the magnetization \M\. The scaling variable \h — h c \L}l v on 
the horizontal axis is the scaled dist ance to the critical point 
while the vertical axis variable, ]M~\ iins L~ d/2 (h - h c ) dv/2 ~ f} is 
the magnitude of the fluctuation in the magnetization, nor- 
malized by the expected number of correlation volumes and 
the magnetization of the correlation volumes of linear size 
£ ~ {h — h c )~ v . The approach to a constant value at large 
scaling variable is consistent with independently oriented re- 
gions for h > h c (upper branch) and magnetization fluctua- 
tions about the ferromagnetic state over regions of size £ when 
h < h c (lower branch.) The solid line represents a power law 
exponent — f3 + dv/2 ~ 1.48, while the upper (lower) dashed 
line has slope 1.63 (1.33.) 

IV. HEAT CAPACITY 

The specific heat of the 3D RFIM is a Quantity 
that can . be measured experimentally, directlyc3E3 or 
indirectly£jc-3 The divergence of the specific heat has 
also been estimated numerically, though not all estimates 
agree and the experimental situation is unclear. Because 
of these discrepancies, it is useful to also study this quan- 
tity in the case of four dimensions, to check the validity 
of the standard scaling picture. 

The heat capacity can be estimated using ground state 
calculations and applying thermodynamic relations em- 
ployed by Hartmann and Young .113 This approach was 
also applied in Refs. [if] and The method relies on 
studying the singularities in the bond energy density 

Ej = L- d Y,s lSj . (6) 

m 

This bond energy density is the first derivative dE/dJ of 
the ground state energy with respect to h (equivalently, 
up to constants, with respect to J.) The derivative of 
the sample averaged quantity Ej with respect to h then 
gives the second derivative with respect to h of the total 
energy and thus the sample-averaged heat capacity C. 
The singularities in C can also be studied by computing 
the singular part of Ej, as Ej is just the integral of C 
with respect to h. The finite-size scaling for the singular 



part of the specific heat C s is 

C s ~L a ^C[(h-h c )L 1 / 1 '}, (7) 

while the scaling for the leading part (through the first 
singular term) of the sample averaged bond energy at 
h = h c is 

Ej, s (L,h = h c ) = Cl +c 2 L^/ v , (8) 

with ci and C2 constants. 

The data analysis is based upon direct fits using Ej. 
This approach avoids complications that arise in com- 
puting the uncertainties when fitting to finite-differenced 
estimates for C, but is otherwise equivalent to fitting to 
such finite differences. The fit was a least squares fit of 
a cubic to h(Ej) for fixed values of L. This fit to the 
inverse function was more stable than fitting to Ej(h). 
The fit function was then inverted to give the estimate for 
Ej(h). The maximum slope of this estimated function is 
in turn used to estimate the peak in C(h) for each L. The 
uncertainties at any point, especially when determining 
the peak value C max (L), in the analysis can be estimated 
using a bootstrap technique (resampling the data.) The 
data for Ej are plotted directly in part (a) of Fig. ||. 
The samples used were the same as used for the magne- 
tization and Binder cumulant analysis. The derivatives 
of the fit are plotted in part (b) of this figure and com- 
pared with the heat capacity values determined by finite 
differencing. Note that the finite differenced values are 
relatively noisy due to the differentiation. This apparent 
noise can be reduced by less refinement in the values of 
h sampled, but this would reduce the resolution in C(h) 
and the location of the peak in C . By directly fitting to 
Ej rather than the finite differences, this complication is 
reduced (but could be managed with appropriate care in 
the error analysis.) 

The estimates for the maximum values of the heat 
capacity are plotted as a function of L in Fig. 0(a). 
The relatively precise data are not consistently fit by a 
power law until L > 16. The fit for these values gives 
ujv = 0.31 ± 0.04, where the errors are purely statis- 
tical. Given the short range of the fit (from L = 22 
to L = 64), one must allow for the possibility of cor- 
rections to scaling giving a different value at larger sys- 
tem sizes (possibly slightly lower.) A fit of these data 
to C max ~ ln(L) is less successful, however (see the in- 
set in Fig. 0(a).) As always, it is difficult to distinguish 
a logarithmic behavior, suggested for C max in Ref. |2^, 
from a small-power-law behavior. In Fig. 0(b), the local 
discrete derivatives are plotted for the cases where the 
behavior should be a power law (main part of the figure) 
and logarithmic (inset.) The power law does seem to be 
more consistent with a convergence to a fixed slope for 
this range of sizes. Though the fits are not definitive, the 
fitted power law behavior is most consistent with scaling 
relations and other data and does seem to explain the 
computed singularity in the bond part of the energy. 

The specific heat can also be used to infer v. This 
can be done directly through scaling the widths of the 
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FIG. 6: (a) The points show the computed dependence of Ej, 
the number density of broken bonds, on h, for L — 4 ... 64 
(not all L values are included, for clarity.) Fits to cubics for 
the inverse function h(Ej) are shown, (b) Estimated heat 
capacity dEj/dh, derived from the differences between the 
E.j values. Solid lines show the derivatives of the fits to Ej 
defined in the text. These derivatives are used to estimate 
the heights of the peaks in the specific heat, C ma x(£). 



peaks in C, but a more robust procedure was to use the 
indirect procedure of fitting Ej, which, being the integral 
over h of C, incorporates the width of the peak in C. The 
quantity (a — was found by using the fitted value of 
Ej at h c . The derivative of Ej(h c ) with respect to ln(L) 
gives the power law (a — = —0.94 ± 0.06. With the 
value for a/v, this gives the estimate v = 0.80 ± 0.06. 



V. STIFFNESS & DOMAIN WALLS 

The nature of responses to external perturbations is 
used to characterize distinct phases, in general. One of 
the more important responses to study is the response to 
changes in the boundary conditions. For example, ferro- 
magnetic phases in pure materials can be identified due 
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FIG. 7: (a) Plot of C max (L) vs. L. The solid line is a fit 
to the finite size scaling form C ma x(£) ~ L~ a ^ , with a/v = 
0.31 ± 0.04. The inset shows a semi-log plot with a least- 
squares fit to the last 3 data points (L > 32.) (b) Local 
derivatives of the plots in (a). 



to the finite energy density of domain walls induced by 
twisted boundary conditions. The application of twisted 
boundary conditions to stiffness and domain walls to 
disordered-jSystems was introduced for spin glasses by 
McMillancj and Bray and Moore. E3 Stiffness and domain 
walls were studied for the 3D RFIM in Rcf. [l(| The ap- 
proach taken quantities studied here for the 4D RFIM 
are the same, though the results are somewhat distinct 
in flavor from the 3D results. 

Measuring the stiffness quantifies the change in en- 
ergy due to a change in boundary conditions. The sym- 
metrized stiffness is defined as 



+ E^ 



E, 



(9) 



where E a b is the ground state energy for boundary spins 
fixed to be a at one end of the sample and b at the other 
end of the sample (periodic boundary conditions are used 
in the other d— 1 dimensions). This definition minimizes 
the effects of surface terms and has the value S = if the 
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two ends of the sample are "decoupled" , with the effect 
of the boundary conditions penetrating only a finite dis- 
tance into the sample. The value E will be zero with high 
probability in the paramagnetic phase, for large samples, 
and is expected to scale as L d ^ 1 in the ferromagnetic 
phase, for fixed h. 



A. Stiffness at criticality 

The sample averaged stiffness E is a quantity that is 
useful for investigating scaling and the order of the transi- 
tion. Near a second order transition, the average stiffness 
scales with a characteristic scale L e , where 9 is the "vio- 
lation of hyperscaling" or stiffness exponent. The natural 
scaling assumption is that this stiffness varies over a scale 
given by the reduced disorder, giving 



Y,nCL 6 S[L x/v (h-hc)K], 



(10) 



with C and K nonuniversal constants and iS a function 
dependent on the shape of the sample. Another charac- 
terization of the distribution of stiffness over samples is 
Po(h,L), which is the probability that the stiffness will 
be exactly zero. As the distribution of the stiffness can 
be scaled at the critical point, with E = invariant under 
rescaling of E by L , Po(h c , L) approaches a constant as 
L — » oo, with the asymptotic value set by sample shape, 
disorder distribution, and lattice type. This convergence 
to a constant was used in Ref. |l6| to locate h c for the 3D 
RFIM. 

The probability of zero stiffness P is plotted in Fig. [|, 
for samples of shape 3L x L 3 . Less anisotropic samples 
had a very small value of Pq and therefore had more 
statistical error. As the running time for a given L is 
larger and the ground states for four different boundary 
conditions were computed, fewer samples were studied 
here than in the magnetization and energy study. For 
L = 32, up to 5 x 10 3 realizations were studied, while for 
the smallest samples, Ej was calculated for about 5 x 10 4 
samples. The estimates plotted are consistent with Pq 
approximately constant in L for h c w 4.18. This is in 
accord with other estimates of h c , though the uncertainty 
in using this plot to determine h c is somewhat larger than 
from other methods. 

A scaling plot showing the collapse of the stiffness cal- 
culations for samples of the same shape, 3L x L 3 , is shown 
in Fig. ^. Assuming the scaling form Eq. ([l0|), a collapse 
to a single function should be found when plotting EL -6 * 
as a function of [h — h c )L 1 l v . This collapse is unrea- 
sonably good (that is, is not too bad for L — 6), using 
h c = 4.177, 6 = 1.82 and v = 0.80. The computations 
strongly support the picture of a secotui.Qrjder transition 
with a value for 9 obeying the boundsuOcJ 



d/2-p/v < 9 < d/2. 



(11) 



When d = 3, it has not been possible to determine 
whether 9 ^ d/2, given that j3/v is so small. Here, given 
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FIG. 8: Plot of the probability of zero stiffness E. The 
samples have a cross section volume of L 3 with a distance of 
3L between the controlled faces. The probability is constant 
to within numerical errors for h w 4.18. 
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the scaled disorder (h — h c )L 1 ^ 1 ' for the values 9 = 1.82 and 
v = 0.80. 



the larger value of 13/ v, it is possible to discriminate be- 
tween 9 and d/2, with the result suggesting that 9 < d/2 
(in addition, the result is consistent with saturation of 
the lower bound.) 



B. Domain walls 

The calculations for the set of boundary conditions 

++, , H — , and — h on the two opposite controlled 

faces (separated by 3L) have also been used to study 
the domain walls in the 4D RFIM. Following the defi- 
nitions of Ref. [l6|, three definitions of the domain wall 
are considered. The first is found by comparing H — (- and 
H — boundary conditions, with the spins fixed to be + 
on the left end in both cases. The set of spins which is 



8 



connected to the left end and fixed under both sets of 
boundary conditions has an internal boundary that in- 
tersects W s bonds. Assuming scaling at h c , the surface 
measure scales as 



W s 



L" 



(12) 



defining the domain- wall dimension d s . The next 
domain-wall measure exponent is found by comparing 
the H — configuration with the and H — h configura- 
tions. The number of bonds which are unsatisfied only 
with H — boundary conditions gives a domain wall mea- 
sure W[. Under and ++ boundary conditions, there 

are unsatisfied bonds due to frozen spin regions, where 
the random field is strong enough to fix the spins under 
all boundary conditions. The unsatisfied bonds with ei- 
ther of these two boundary conditions are not counted as 
part of the H — domain wall under this definition. The 
only broken bonds which are counted as part of Wj are 
those broken due to the twisted boundary conditions. 
This measure similarly defines an incommensurate sur- 
face exponent by 



Wj 



(13) 



The third definition of the effect of boundary conditions 
is given by the bond or exchange part of the stiffness, 



£ 7 = (K- 



E J _ + -E{ + -Et_)/2, 



(14) 



J^2(ij\ SiSj. This count includes some bro- 



where Ej 

ken bonds with negative sign and is influenced by frozen 
islands. The "dimension" dj is then 



L" 



(15) 



As Ej is the deriwitive of £ with respect to J, thermo- 
dynamic relational!] imply that 



e + l - 

V 



(16) 



The values of d s , dj and dj were estimated by taking 
the discrete logarithmic derivatives of W s , W / and Ej, 

d s ,i,j(VLiL 2 ) = In[W(L 2 )/W(ii)]/ ln(L 2 /Li), (17) 

with W being one of the measures of the domain wall. 
The results are plotted in Fig. [H]. 

One of the more striking differences between the 3D 
and 4D calculations is that, while d s ^ d in 3D, when 
d = 4 the value of eL is consistent with the relation 



gL = d. 



(18) 



Thus, the domain wall defined by the surface of connected 
fixed spins anchored at the fixed end of the sample has 
dimension consistent with the spatial dimension. This 
surface, the internal boundary between flipped and fixed 
spins, appears to be space-filling. Additionally, the esti- 
mated value of di = 3.20±0.12 is clearly distinct from d s , 
in contrast with the near equality seen in 3D.E3 These dif- 
ferences will be addressed in more detail in Sec. VI, when 
examining the frozen spins for h < h c . 
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FIG. 10: Estimate of the dimensions (a) d a , (b) di, and (c) 
dj obtained from the discrete logarithmic derivatives of wall 
area W a , the number of bonds Wi created by twisted BC's 
relative to uniform sign BC's, and the exchange stiffness Ej 
with respect to L for several h. These plots are used to infer 
d a = 3.94 ± 0.06, di = 3.20 ± 0.12 and dj = 2.94 ± 0.12. 



C. Domain walls and scaling 

The value computed here for dj, dj = 2.94 ± 0.12, 
is just consistent with Eq. Jl6|). As the derivation of 
Eq. ( |l6| ) is quite robust, this consistency should not be 
surprising. Though the arguments arerapparently sound, 
the conjecture made originally for 3Dji3 namely that 



d,, + 0/v, 



(19) 



9 



has a non-rigorous derivation, especially as in the above 
form, d s in the original version has been replaced by 
di, which here more clearly reflects the measure of the 
domain walls induced by boundary condition changes. 
However, this scaling relation is easily consistent with 
the computed values of the domain wall exponents and 
P/v. 

The value found here for (a — l)/v also satisfies the 
relationship dj — d — f3/v = (a — which was used 

in Ref. |l6|, except for the replacement here of d s by dj, 
motivated in the 4D RFIM by the more natural defini- 
tion of domain walls using di and the spatial structure 
of frozen spins. 



VI. FROZEN AND MINORITY SPINS 

The result that one measure of the domain wall dimen- 
sion, d s , is near to the spatial dimension d suggests that 
the picture of the spin configurations must differ between 
the cases d = 3 and d = 4. The picture of the configu- 
ration at the transition in d = 3 described in Ref. [lq is 
that of nested domain walls, where the domain walls are 
the boundaries separating connected sets of spins of the 
same sign (see also Ref. [19].) In this section, results are 
presented that necessitate a different picture in d = 4, 
due to the percolation of minority spins for h less than 
the critical disorder h c . (These results should be com- 
pared with those for the 3D RFIM presented in Ref. |28| , 
which support the existence of two interpenetrating span- 
ning domains in the 3D RFIM for h > h c , and those of 
Ref. ^9|, where the surprising claim is made that there is 
a second critical h p > h c where there is first simultaneous 
spanning by up and down spin clusters.) The percolation 
of minority spins in d = 4 for h < h c makes the identi- 
fication of domain walls with connected sets of uniform 
spins problematic. 

Given a disorder realization {hi}, there are two natural 
sets of spins to consider when defining domain walls and 
percolation clusters. The minority spins are simply those 
that have spin opposite to the mean magnetization. The 
fraction of spins that fall into this category is (1— |m|)/2. 
Frozen spins are those that are invariant under all bound- 
ary conditions. These spins are minority spins under ei- 
ther all up or all down boundary conditions, so that the 
fraction of frozen spins is 1 — \m\, when L>^. If either 
minority or frozen spins were distributed independently 
in space, the clustering of these spins would map directly 
onto simple percolation. As there are strong interactions 
between these spins and the boundaries between them 
are related to domain walls, the percolation is not sim- 
ple, on short length scales. For h < h c , the correlations 
should vanish in the limit of separations much greater 
than £. 

The clustering and percolation behavior of these spins 
can be directly studied to learn more about the do- 
main walls. From ground state configurations for L = 8 
through L — 64, computed both for all up and all down 



boundary spins, the frozen and minority spin sets were 
identified. These sets can be studied directly or in a 
coarse grained sense. (Coarse grained spin blocks were 
determined by whether the minority or frozen spins were 
a majority of the block, with ties randomly broken.) The 
spanning clusters were defined as those that connected 
two opposite faces of the hypercubic sample. 

Fig. is a plot of the percolation probability p b (i.e., 
the probability of at least one spanning cluster) of minor- 
ity spins on scales b — 1, 2, 4 as a function of h. From this 
plot, an extrapolation of the curve crossings to large L 
suggests that the minority spins percolate in the infinite- 
volume limit at h™ = 3.850 ± 0.005. Plots for the frozen 
spins are qualitatively similar, with a lower percolation 
threshold of h{ = 3.680± 0.005. In each case, the number 
of spanning clusters peaks near h™'*, with the peak num- 
ber increasing with L. The percolation point tends to- 
ward h c as the scale b increases, consistent with a scaling 
toward the ferromagnetic state of uniform magnetization 
for h < h c . 

Implications for the simple domain wall picture in 
d = 4 follow directly. As the minority spins percolate 
in the ferromagnetic phase, h < h C) the boundaries of 
connected sets of same-sign spins are space filling. The 
definition of domain walls using simply these connected 
sets is thus not clearly informative about the effect of 
boundary condition changes. The surfaces defined by the 
incongruent bonds are more useful in understanding the 
domain walls. These are the surfaces that separate the 
two ferromagnetic ground states; the frozen spins make a 
space-filling background that is common to both states, 
even for h < h c . (In three dimensions, at h — /i c , there is 
a fractal set of spins that can be controlled by the bound- 
ary conditions.) In addition, the relationship in 3D be- 
tween (3/v and poo, the probability of crossing a domain 
wall per factor of e in length scale, would be much more 
difficult to investigate in 4D, as the domain walls are not 
readily identifiable. 



VII. STATES 

In earlier sections, it has been (sometimes implicitly) 
assumed that the transition in the 4D RFIM is consis- 
tent with the simple picture of a ferromagnetic to para- 
magnetic transition, with the sign of the magnetization 
in the ferromagnetic state dependent on boundary condi- 
tions and the spin configuration independent of boundary 
condition in the paramagnetic state, far from the bound- 
aries. This assumption is examined in this section. .The, 
approach is inspired in large part by analytic work.EHEU 
A discussion of the numerical study of the nature of the 
thermodynamic states is presented in Ref. ^l] and the ap- 
plications to the 3D RFIM can be found in Ref. In 
summary, one test of the number of states in the thermo- 
dynamic limit is to determine the correlation functions 
(in this case, the ground state) in the interior of the sam- 
ple under several different boundary conditions. For a 
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FIG. 11: Plot of the percolation probability for minority spins 
as a function of disorder strength h, for L — 8, 16, 32, 64. The 
dashed line indicates h = h c . (a) Percolation probability pi 
for minority spins. The percolation probability approaches 
unity for large systems for hp — 3.875 ± 0.005. (b) Percola- 
tion probability p2 for minority block spins of size 2 4 . The 
percolation threshold is closer to h c . (c) Percolation proba- 
bility p4 for minority block spins of size 4 4 . 



uration in the interior of the sample. In particular, the 
periodic (P), all spins up (+), all spins down (— ), and 
open (O) boundary conditions were compared. As the 
4D computations are much more time consuming, the 
comparison between the ground states of a system and a 
smaller subsystem, each with open boundary conditions, 
was not extensively studied, as it was in the 3D RFIM. 

A summary of the results for comparisons between P 
and +/— is presented in Figs, [l^ and 13. The results for 
O vs. +/— are quite similar. The plots show the (scaled) 

probability Pp.^ (2, h, L), for a given h and L, that the 

ground-state P configuration is distinct from both the 
ground-state + and — configurations in the central vol- 
ume of size 2 4 . As L increases, this probability decreases 
toward zero at all h. This suggests that in large sam- 
ples, the interior configuration for a number of boundary 
conditions (including periodic and open) can be found by 
imposing either + or — boundary conditions. It was also 
found, as in the 3D case, that for h > h c , as L increases, 
the interior configurations for the + and — boundary con- 
ditions become identical with unit probability. Together, 
assuming the extrapolation to large L is correct, these 
results show the existence of a single state for h > h c 
(for if the interior configuration differs between any two 
boundary conditions, it must differ between + and — ) 
and strongly suggest the existence of only two states for 
h < h c . 

The scaling of P P j— is. consistent with previous work 
on disordered modelsOS As a function of the scaled dis- 
order (h—h c )L 1 / u , the function approaches a single curve 
when the probability Pp^ is scaled by L d ~ dl . This scal- 
ing results from assuming that the number of large (size 
L) domain walls induced by generic boundary condition 
changes is constant as L — * oo. This assumption is con- 
sistent with the observation that d] > d— 1. The chance 
that an interior volume of fixed linear size w intersects 
a domain wall is expected to behave as (w/L) d ~ dl . The 
clean collapse of the data for larger system sizes, using 
values determined from magnetization and domain wall 
measurements, lend quantitative support to this picture. 

The dimension di could be alternatively deduced from 
this data. Fig. |lj shows the dependence of the peak value 

of Pp,-\ on L. Assuming d — di gives this slope, the data 

for L = 12 -> 44 gives rfj = 3.19(2). This value could be 
taken as the best one for di, but this is not done here and 
is instead used as a confirmation of the scaling picture. 



small number (one or two for Ising models) of thermo- 
dynamic limits, there will be a small number of interior 
configurations. The probability that the interior of the 
ground states will differ from one of the large volume 
limit configurations decreases as a power law dependent 
on the dimension of the domain wall. 

The degeneracy of the ground state was directly ad- 
dressed for the 4D RFIM by studying the effect of chang- 
ing boundary conditions on the ground state spin config- 



VIII. ALGORITHMIC SLOWING DOWN 

In the 3D RFIM, it was found that the number of oper- 
ations carried out by the ground state algorithm diverged 
near the ferromagnetic-paramagnetic transition. In this 
section, similar results are presented for the 4D RFIM. 
The scaling results here have greater accuracy near the 
transition. The scaling is found to be quite consistent 
with the heuristic picture presented in Ref. |3^. 

The key quantitative relation to be elucidated is that 
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FIG. 12: Plot of P P ,+_(2, /t, L), the probability that, at given 
sample size L and disorder h, the ground state with periodic 
boundary conditions will differ from ground state configura- 
tions for both fixed + and — boundary conditions, in a volume 
of size 2 4 in the center of the volume L 4 . Note that extrap- 
olating to L — > oo suggests that Pp^ (2, h, L) — * for all h. 

The solid lines are least-squares fits to quartics in ln(P). 
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FIG. 13: Scaling plot for Pp^ (2,h,L), using the scaled 

probability Pp i +_L d ~ dj and scaled disorder (h—h^L 1 ^, with 
values di — 3.2 and v — 0.8. The scaling collapse suggests 
that the changes in boundary conditions typically introduce a 
finite number of domain walls (say, one of size L) with ~ L dl 
bonds at h = h c . 



between the number of primitive operations carried out 
to find the ground state and the physical understand- 
ing of the phase transition and correlation volumes. In 
Ref. it was argued that the time per spin to find the 
ground state in the RFIM is directly proportional to the 
linear size L near the transition. Tbisjesults from the 
nature of the push-relabel algorithmtS'Eil used, which ef- 
ficiently constructs a "height" field over the lattice that 
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FIG. 14: Plot of the dependence of the maximum of 

Pp t -I (2, h, L) over h on the system size L. The solid line 

is a fit to PP^_ ~ L- (d - dl \ with d-d! = 0.81 ± 0.02. 
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FIG. 15: (a) Plot of the number of relabel operations per spin, 
r, carried out by the ground state algorithm as a function of h, 
for L = 4 . . . 64. (b) Plot of the sample-to-sample fluctuations 
o r in the number of relabel operations per spin. This quantity 
provides an especially sharp and quickly diverging curve for 
estimating h c . 
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FIG. 16: Linear plot of the height of the peak r p (L) in r vs. 
h. The value of r p scales almost exactly linearly with L. 

guides the "flow" (corresponding to the redistribution of 
excess "fluid" or external magnetic field hi.) 

The relabel operation is one of the primitive operations 
carried out during the convergence of the algorithm to the 
physical ground state. In Fig. ^[ the sample average r of 
the number of relabel operations per spin and the sample- 
to-samplc fluctuations in r, ay, are plotted as a function 
of h for different L. (Results for the number of the other 
primitive operations, the push operations, are quite sim- 
ilar.) There is clearly a peak in both quantities near 
h c , with the peak in the sample-to-sample fluctuations 
being more sharply peaked, relative to the non-critical 
contribution^. 

The r(h) curves were fit with fourth-order polynomials 
to extract the peak value r p (L). The plot of this quan- 
tity is shown in Fig. A linear fit is shown, which is 
remarkably consistent with the data over a wide range 
of L. The result is in agreement with the arguments of 



Ref. |33| and supports the relationship between the physi- 
cal correlation length and the evolution of the algorithm. 



IX. SUMMARY 

By computing the ground state for a large number of 
samples of volume up to L 4 = 64 4 , the quantitative and 
qualitative thermodynamic properties of the 4D RFIM 
have been studied. The derived exponents satisfy the 
conventional scaling relations. The values of the expo- 
nents and location of the transition are consistent with, 
but are based on larger systems than, the results for the 
Gaussian 4D RFIM published in Refs. |ll|,||. Note that, 
as in previous work, the value for v is the least certain 
and that errors in v propagate to estimates of j3 and a. 
The picture of a single transition from a nearly-two-fold- 
degenerate ferromagnetic state to a single paramagnetic 
state is confirmed by comparing ground states with vary- 
ing boundary conditions. Strong evidence-is presented 
that the picture of domain walls developedlla for the 3D 
RFIM must be modified to describe the 4D RFIM. In 
particular, the percolation of frozen and minority spins 
within the ferromagnetic phase implies that the sets of 
connected same-sign spins are not the boundaries of do- 
main walls. The empirical running times for the ground 
state algorithm peak near the phase transitkfflr,in a man- 
ner consistent with previous descriptions,t3E2l with the 
peak running time per spin apparently proportional to 
the linear system size. The algorithmic running times 
provide a check on the location of the transition and the 
scaling exponent v. 
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